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Abstract 

Shear transformations (ie., localised rearrangements of particles resulting in 
the shear deformation of a small region of the sample) are the building blocks 
of mesoscale models for the flow of disordered solids. In order to compute the 
time-dependent response of the solid material to such a shear transformation, 
with a proper account of elastic heterogeneity and shear wave propagation, we 
propose and implement a very simple Finite-Element (FE) - based method. 
Molecular Dynamics (MD) simulations of a binary Lennard-Jones glass are 
used as a benchmark for comparison, and information about the microscopic 
viscosity and the local elastic constants is directly extracted from the MD 
system and used as input in EE. We find very good agreement between EE and 
MD regarding the temporal evolution of the disorder-averaged displacement 
field induced by a shear transformation, which turns out to coincide with the 
response of a uniform elastic medium. However, fluctuations are relatively 
large, and their magnitude is satisfactorily captured by the FE simulations of 
an elastically heterogeneous system. Besides, accounting for elastic anisotropy 
on the mesoscale is not crucial in this respect. 

The proposed method thus paves the way for models of the rheology of 
amorphous solids which are both computationally efiicient and realistic, in 
that structural disorder and inertial effects are accounted for. 

Keywords: shear transformation, plastic event, structural disorder, elastic 
moduli 

PACS: 62.20.D-, 83.80.Ab, 02.70.Dh, 61.43.Bn 


1. Introduction 

Glasses are macroscopically isotropic and homogeneous. Microscopically, 
the absence (or elusiveness) of a clear structural signature of the liquid-to- 
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glass transition upon cooling may fallaciously lead one to believe that the 
property of structural homogeneity holds down to the microscale in glasses, 
as it does in liquids. In the past decades, it has become clear that this idea 
is completely erroneous: structural disorder is often regarded as a crucial 
aspect not only of the glass transition, but also of the flow of soft or hard 
glassy materials, and more generally amorphous solids, such as emulsions, 
foams, dense gels, and granular matter. For instance, several theories of the 


glass transition, including the Random First Order Theory (Lubchenko and 


Wolynes, 2007; Berthier and Biroli, 2011) and kinetically constrained models 


with facilitated dynamics ( Chandler and Garrahanf 2010), put the focus on 
dynamical heterogeneities, that is, the coexistence of regions with fast and 
slow (arrested) dynamics. 

Heterogeneities are even more manifest when the materials are forced to 
flow. Instead of a homogeneous deformation, one observes localised bursts 
of particle rearrangements, called shear transformations or plastic events. 


embedded in an essentially elastically deforming medium (Argon and Kuo 


1979[ [Falk and Langer[ |1998[ [Schall et aLf |2007t |Amon et akf |2012| ) . These 
irreversibly rearranging regions coincide with “weak” zones where the local 


elastic (shear) moduli vanish at the onset of a plastic event (Tsamados et al. 


2009). By simply looking at the instantaneous (static) configuration of the 


system, computing its soft modes, and observing where they concentrate, one 
can predict statistically (but only to a limited extent) the position of future 
rearrangements (Widmer-Cooper et ah, 2008; [Rottler et aT , 2014). Not only 
does microscopic structural disorder play the leading role in fixing where the 
rearrangements will occur, but it also affects the way stress is redistributed 
in the medium during these plastic events, i.e., the propagation of the shear 
waves originating from the rearranging region. 

This stress redistribution is generally described as the solution of an Es- 
helby inclusion problem in a uniform linear elastic medium (Eshelby, 1957), 
with an inclusion that is often assnmed pointwise in lattice-based rheolog- 


ical models, for convenience (Picard et al., 

2004 

, |2005l 

Vandembroucq and 

Roux, 2011; Talamali et ah, 2011| Lin et ah, 

2014[ Martens et al. 

, 2011, 

2011 

Nicolas et ah, 2014a). The solution is given by an elastic propagator Q with 


a characteristic four-fold angular symmetry and an r ^ spatial decay in two 


dimensions, in line with experiments on, e.g., dense emulsions (Desmond and 


Weeks, 2013) (also see Budrikis and Zapperi (2013); Sandfeld et al. (2015) for 


a discussion on this elastic propagator and its possible numerical implementa¬ 
tions). However, some of us very recently showed that such description only 
holds on average ( Puosi et al.| 2014); if an individual plastic event is con¬ 
sidered, the description is unreliable, because the average response is blurred 
by sample-to-sample fluctuations, presnmably associated with the elastic het¬ 
erogeneity of the material. Moreover, this approach neglects inertial effects 
by supposing instantaneous mechanical equilibration, or, in other words, an 
infinite shear wave velocity, whereas the role of inertia on the statistics of 


avalanche sizes has been numerically evidenced (Salerno et ah, 2012; Salerno 


and Robbins, 2013). These two deficiencies, possibly among others, under¬ 


mined a recent endeavour of ours to reproduce the spatio-temporal correla¬ 
tions in the flow of a disordered solid with a coarse-grained model using the 
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elastic propagator Q (Nicolas et al., 2014b). 

The objective of this contribution is to go beyond the average, equilibrium- 
based description in terms of the elastic propagator; we aim to devise and 
put to the test a minimal framework allowing to capture the fluctuations in 
the response due to structural disorder, as well as the propagation of the 
shear waves, in two dimensions (2D). To this end, we implement a basic 
Finite Element (FE) code and use Molecular Dynamics (MD) simulations of 
an athermal solid as a benchmark. In so doing, we show how the microscopic 
data about, e.g., the local elastic constants can be extracted from the MD 
system and used as input in FE. 

In Section 2, we present the MD simulation method and we introduce our 
simplified FE algorithm. Section 3 is concerned with the fitting of the me¬ 
chanical parameters required by FE, in particular, the calculation of the local 
elastic constants of the MD solid. Section 4 clarifies the protocol to trigger 
artificial shear transformations. Finally, Sections 5, 6, and 7 describe the 
disorder-averaged elastic response to this localised transformation, the fluc¬ 
tuations around this average, and the response in a particular configuration 
of the system, respectively. 


2. Methods 

2.1. Molecular Dynamics 

To probe the fiow properties of amorphous solids, we resort to MD sim¬ 
ulations of a 2D amorphous system. More precisely, we simulate a binary 
mixture of A and B particles, with = 32500 and Nb = 17500, of re¬ 
spective diameters gaa = 1-0 and ass = 0.88, confined in a square box of 
dimensions 205a aa x 205(J^^, with periodic boundary conditions. The sys¬ 
tem is at reduced density 1.2. The particles, of mass m = 1, interact via a 
pairwise Lennard-Jones potential, 


14/3 (r) = 4e„/3 




5 


where a,/3 = A, B, gab = 0.8,e^^ = 1.0, cab = 1-5, and cbb = 0.5. The 
potential is truncated at r,. = 2.5cr^^ and shifted for continuity. 

We conduct our study in the athermal limit, by thermostatting the system 
to zero temperature, so that no fluctuating force appears in the equations of 
motion, viz., 


m- 


dvj 

dt 

dvj 

dt 


= Vi 




dV (ry) 


dr. 


+ 


D 


^3 


( 1 ) 


The dissipative force fi^ experienced by particle i is computed with a Dis¬ 
sipative Particle Dynamics (DPD) scheme, whereby particles are damped on 
the basis of their relative velocities with respect to their neighbours. More 
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precisely, fi^ reads 


fi^ = 
where w(r) = 


inj) 

j^i ij 

1 - ^ if r < Tc, 

0 otherwise. 


( 2 ) 


Here, Vij = Vi — Vj denotes the relative velocity of particle i with respect 
to j, the vector Vij = r* — Vj connects particle j to i, the cut-off distance is 
set to Tc = 2.5aAA, and ^ controls the damping intensity. Different values of 
( will be tested to probe the different damping regimes, from underdamped 
(C ^ 1) to highly overdamped (C ;:§> 1). Note that, in Eq. the projection 
of the force onto the radial vector Vij is required in order to conserve angular 


momentum. Several other virtues of DPD have been exposed by |Soddemann| 
et ah (2003). As far as we are concerned, one of the main advantages is that, 


in the light of the recent work of Varnik et al. (2014), experimentally measured 


correlations in the flow of amorphous solids are better reproduced numerically 
when dissipation is based on relative particle velocities, in opposition to a 
mean-field damping scheme, in which absolute velocities (with respect to a 
hypothetic solvent flow) are used. The impact of this implementation on the 
propagation of shear waves will be discussed in Section |5.2[ 

However, the DPD algorithm does not conserve the position of the centre 
of mass of the system a priori. Since the ensuing global translations of the 
system may disturb the forthcoming analysis of displacements in reponse to 
shear transformations, the system is regularly re-centred during the simula¬ 
tion. 

Equations [T] are integrated with the velocity Verlet algorithm with 5t = 


0.005. In all the following, we use tlj 
uaa as the unit of length. 


\lma\Ale as the unit of time and 


2.2. Simplified Finite Elements 

In the presence of elastic heterogeneities, the elastic response to a lo¬ 
calised shear transformation becomes intractable to analytical calculations. 
This notably implies that the Fast Fourier Transform routine commonly used 
in elastoplastic models needs to be replaced. As a minimal substitute, we 
propose a simplified FE algorithm, which will also allow us to account for 
inertial effects. 

The FE method consists in discretising a Continuum Mechanics equation 
onto a mesh. Here, the Continuum Mechanics equation involves elastic and 
dissipative (viscous) forces, as well as inertia; hence, the momentum conser¬ 
vation equation reads 


p—(T’,t) = y • [C{rfi)e{rfi)] + ? 7 V^A(r,t) , 

^ r ^ elasticity viscosity 

inertial force 


(3) 


where u and e are the displacement and strain fields, respectively, ^•/ot = 
^•/dt+{v ■ V) • denotes the convected derivative, dots denote time derivatives, 
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Figure 1: Sketch of the FE mesh. The system is periodic in both directions, 
so that column coincides with column 0 and row Ny coincides with row 
0. There are N = x Ny nodes and elements. 


p is the (area) density of the material, C denotes a local stiffness matrix, and 
r} is the microscopic viscosity. Upon discretisation, it turns into 


• il = 

inertial force elasticity viscosity 

where u is now a shorthand for the high-dimensional vector 

T 


(N-1) 
U)c h 


U 


(N-1) 

y 5 


( 0 ) 
U)c h 


U 


( 0 ) 


(4) 


containing the displacements along x and y at the N nodes of the mesh. 
A4, /C, and H are 2N x 2N real matrices (to be specified later), and the 
dependences on time have been omitted. 

Bearing in mind our pursuit of minimalism, we choose a simple (static) 
regular square meshgrid, as sketched in Fig. [T} In an element, the local strain 
€ = ( €xx, £yy, V^^xy ) ) using condcnscd notations for 2D symmetric ten¬ 
sors, is a function of the displacements at the local nodes, and we make the 
approximation of a uniform strain within each elemenlj^ For convenience, let 
us number these nodes from 0 to 3 counter-clockwise, for a given element, 
starting from the bottom left corner, viz., In an analogous way, the 

(uniform) elemental stress cr®* is derived from the nodal forces (/j', /®'). Since 
the mesh is regular, we can dehne a constant 3x8 real matrix B that relates, 
in a given element, the (nodal) displacements to the (elemental) strains, on 
the one hand, and the (nodal) forces to the (elemental) stresses, on the other 
hand, viz., 


^ In practice, our simplified FE method is therefore close to a Finite Volume method. 
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The expression of the matrix B is given in [Appendix A , along with further 
details pertaining to the implementation of the FE routine and the compu¬ 
tation of the matrices A4, /C, and H appearing in Eq. Note that the 
\/2 prefactors have been introd uced with foresight (see Section | 3 . 2 | ) and the 
“minus” sign preceding B in Eq. A.l is due to the fact that ^ is the force 
exerted by the element on node i. 

The resulting routine is still simple enough to be used quite efficiently in 
a coarse-grained model. In particular, (see Appendix A.4), the global force- 


displacement matrix is constant and, accordingly, only has to be inverted 
once, at the beginning of the simulation. 


On the other hand, there are naturally a few downsides to this simplicity. 
Eirst and foremost, it is only marginally stable, insofar as the convergence 
of the discrete FE solution to the continuous solution of Eq. is not guar¬ 
anteed when the mesh size tends to zero. Consequently, this scheme is not 
suited to general purpose. However, as will be shown below, it is both sat¬ 
isfactory and very convenient for the modelling of (the response to) shear 
transformations, where elements represent material regions of hnite size. In 
particular, the frequently encountered checkerboard issue, whereby high and 
low displacements/velocities alternate erratically in neighbouring cells (hence 
the image of a checkerboard), is practically circumvented, provided that shear 
transformations span four adjacent elements (a “macro-element”) and inertia 
is present, i.e., p ^ 0. 


3. Fitting of elastic and viscous parameters 

We are now left with the task of fitting the physical parameters appearing 
in Eq.|^with the MD parameters. Neglecting mesoscopic density fluctuations, 
the density p and the miscroscopic viscosity p are supposed to be constant, 
while the stiffness matrix C{r,t) is allowed to vary in space. 


3.1. Viscosity 

To fit the viscosity p in Eq. we compare the stress due to homogeneous 
shear, at a rate 7 , as calculated, on the one hand, in FE (a^y = P'j), and, 
on the other hand, in MD (where it is obtained through the Irving-Kirkwood 
formula). The calculations are shown in their full extent in Appendix B and 


lead to the following formula for a binary mixture of A and B components: 


r 

3= -X 


POO 

/ [n\gAA{r) + 2nAnBgAB{r) + n%gBB{r)] (r) r^dr, 
Jo 


7 












where ha and ns are the number densities of A and B constituents in the 
system, qaa, Qbbi and qab are the radial distribution functions for the A — A, 
B — B, and A — B correlations, respectively, and ^ and w are the DPD 
damping coefficient and the damping function dehned in Eq. 

For the MD system under consideration, we obtain 

7] = 0.726 C- 


3.2. Local elastic constants 

Having determined the dissipative coefficient of the model, we turn our 
attention to the local elastic properties of the system. 

The only relevant material lengthscale in the model being the typical 


size (a = 5a aa) of a rearrangement (Nicolas et ah, 2014b), we tile the sys¬ 


tem into subregions of size a and compute the local stiffness tensors on this 
“mesoscopic” scale, with the local stress-affine strain method presented in 


Ref. (Mizuno et ah, 2013). Details of this protocol and issues related to the 
rather unfamiliar local stiffness tensors are discussed in [Appendix C[ With 


condensed notations, these tensors can be written as 3 x 3 real matrices in 
2D, viz., 


a. 


'yy 


\/2c 


xy 


a 

= I c, 
a 


yy,xx 


xy,xx 


a 

a 

a 


XX,yy 

yy^yy 


xy,yy 


a 

a 

a 


XX,xy 

yyay 


^yy 


xy,xy 


y/2t 


( 6 ) 


xy 


where a^ 


Gyy, and a^y are the linear elastic contributions to the local stress. 


Contrary to their macroscopic counterpart, the local C matrices are not 


symmetric a priori, for very small regions (Tsamados et ah, 2009). However, 


the coarse grain a = 5a aa is large enough here for the assumption of symme¬ 
try to be a reasonable approximation. To limit the number of parameters, we 
further assume that isotropic contraction/dilation of the region only generates 
an isotropic stress, i.e., that 


( 


^yy 


\/2e, 


xy 


= C 2 / 2 ( 1 1 0 ) 


T 


is an eigenvector of C. 

These two assumptions, namely, tensorial symmetry and isotropy of the 
response to contraction, imply that the stiffness tensor should be of the form 


/ a 5 \ 

C =[ 5 a -(3 \, (7) 

V/5 -/5 V J 


where the parameters a, 5, /?, u G M are assessed in [Appendix C[ By analogy 
with the macroscopic situation, the eigenvalues ci ^ C 2 ^ C 3 of the approxi¬ 
mated matrix C are related to the local shear moduli ni and 112 and the local 
bulk modulus K via ci = 2/ri, C 2 = 2 /^ 2 , and C 3 = 2K, and there exists a 
frame {e^iO), ey{6)), rotated by an angle 0 with respect to the original frame, 
in which the stiffness tensor reads 


K -\- [12 LC — P 2 0 

K — 1^2 K -\-1^2 0 

0 0 2 pi 


with /ii ^ /i 2 - 



















Denomination 

Symbol 

Mean 

Std dev. 

Shear modulus (weak direction) 

jii 

13.16 

7.2 

Shear modulus (strong direction) 

S2 

24.46 

5.8 

Average shear modulus 

— M1+M2 

— 2 

18.81 

5.3 

Bulk modulus 

K 

99.9 

8.4 


Table 1 : Statistical properties of the elastic constant distributions: mean 
values and standard deviations (std dev.). 


Consequently, the following four local parameters suffice to determine C com¬ 
pletely: 6, pi, fj, 2 , and K. 

Table [T] summarises the main features of the distributions of pi, /i 2 , and 
K measured in the Lennard-Jones glass under consideration; 9 is uniformly 
distributed, in accordance with macroscopic isotropy. 

It is noteworthy that the local stiffness matrices exhibit significant anisotropy, 
as indicated by the discrepancy between the mean value of the shear modulus 
in the (locally) weaker direction, (pi) = 13.16, and its strong counterpart, 
(/i 2 ) = 24.46. 

Some regions actually even display negative shear moduli pi. This is 
not unrealistic in the MD system, because these regions can be stabilised by 
the surrounding medium, but in the following they will be discarded, and 
arbitrarily set to zero, in the FE simulations, where they cause instabilities. 

Lastly, the bulk modulus is much larger (by a factor of 5) than the shear 
moduli, in line with expectations, and its relative standard deviation (he., 
the ratio of the standard deviation and the mean value) is by far smaller 
than it is for the shear moduli, which means that, on a relative basis, the 
latter are more broadly distributed. Consequently, we will henceforth always 
neglect spatial fluctuations of the bulk modulus and set iF = 99.9. As for 
the distributions of shear moduli, three types of systems will be considered 
in FE: 

(i) a uniform system, with = P 2 = 18.8 

(ii) a heterogeneous system made of isotropic blocks (“het. iso.”), with 
1^1 = 1^2 = 18.8 ±5.3, he., a normal distribution of shear moduli pi = /i 2 with 
mean value 18.8 and standard deviation 5.3. (Remember that each block is 
a macro-element made of four adjacent finite elements.) 

(in) a heterogeneous system made of anisotropic blocks (“het. aniso.”), 
with /ii = 13.16 ± 7.2 and /i 2 = 24.46 ± 5.8 and a uniform distribution of the 
angles 9. 

Through the simulation of plane shear waves, we have checked that the 
transverse sound velocity measured in FE is consistent with that measured 
in MD. 

4. Protocol for the artificially triggered shear transformations 

In this section, we describe the protocol to artificially trigger ideal shear 
transformations. 
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Figure 2: Sketch of the displacements applied to a macro-element to model a 
pure shear transformation. 


In the MD system, following Puosi et al. (2014), shear transformations 
are artificially created by applying a pure shear strain to a disk centred at 
and of diameter a = 5cr^^. To do so, particles whose initial position 
{xi,yi) belongs to this region are moved to a new position at t = 0, 

which satisfies 


j Xi —>■ x'^ — Xi + €xy {yi — yo) 

\yi ^ y'i = yi + ^xy ~ ^o) • 


Their positions are then frozen for the whole simulation. In order to mea¬ 
sure the elastic, he., reversible, response of the medium, e^y never exceeds 
a few percent strain. Clearly, all (transient or permanent) dilational effects 
(Schuh et ah, 2007) potentially accompanying shear transformations are here 
discarded. 

A similar shear transformation is applied in the FE simulations to a macro¬ 
element made of four adjacent elements (see Section 2.2), by controlling the 
positions of the nodes of these elements, as sketched in Fig. 


5. Disorder-averaged propagation of shear waves 

Let us first probe the disorder-averaged time-dependent response to a 
shear transformation, in different damping regimes, both in FE and in MD. 
To this end, MD simulations are averaged over many (50) locations of the 
shear transformation in the sample, while the FE results are averaged over 
many (50) realisations of the disorder, i.e., of the random values of the local 
elastic constants. 


5.1. Comparison between MD and Finite Elements 

Eor a quantitative study, we make use of the average propagation radius 


Ar(t) introduced by Puosi et al. (2014) to measure the advance of the wave. 


^r{t) = // \ur{rC)\d^r, 


where Ur{t) is the radial displacement at time t. If the final displacement 
{ur{r]t = oo) ~ r~^ in any given direction 0 in the far field) is essentially 
achieved as soon as a region is reached by the wavefront, Ar{t) will grow 
linearly with the (linear) size of the displaced region. The average propagation 
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radius is plotted in Fig. for diverse values of the damping The initial 
growth is ballistic in MD, with Ar{t) ~ t , while at long times Ar(t) saturates 
to its steady-state value. The evolution of Ar{t) before the steady state is 
reached strongly depends on (^. At low damping (C = 1), the interaction 
with the waves generated by the periodic replicas of the shear transformation 
leads to particularly long-lived oscillations of Ar{t) (Fig. [^), while stronger 
damping = 100) completely suppresses these oscillations. 

The FE simulations nicely capture this qualitative change, and the agree¬ 
ment both in the limit of low damping (Fig. [^) and in the limit of strong 
damping (Fig. [^) is excellent, at relatively long times. This is true for all 
three FE systems, including the uniform one, which snpports the idea that the 
average propagation in elastically heterogeneous media is virtually identical 
to the propagation in a uniform medium. 

Eor an intermediate value of the damping, namely, ^ = 10 (Eig. [^), 
the agreement is reasonable, but not quite as good, insofar as the oscillations 
observed in MD are damped perceptibly faster than their connterparts in FE, 
not only in the uniform system, but also in the heterogeneous one (het. iso.). 
This suggests that the EE viscosity is somewhat underestimated, or that the 
anharmonicities present in MD significantly contribute to the damping of the 
oscillations. 

Einally, the short-time propagation is well described at low damping, but 
the agreement declines when C increases, in which case the FE method over¬ 
estimates the propagation velocity over short distances. 


5.2. Theoretical rationalisation 

Puosi and co-workers (Puosi et al. , |2014 ) reported that, with a mean-field 
dissipative force (i.e., by substituting fi = for Eq. [^in Eq. [^, Ar{t) 

initially grows in a diffusive fashion, i.e., Ar{t) ~ at large damping, 
that is to say, for short Langevin damping times r^. By contrast, no such 
diffusive regime is observed here, even for large damping parameters (. The 
dissipation scheme therefore affects the nature of shear wave propagation. 
Can this discrepancy be explained theoretically? 


5.2.1. Mean-field dissipation 

In the presence of a mean-field damping force, force balance on particle i 
can schematically be written as 

mvi{t) ^ iuj{t) - Ui(t )), ( 8 ) 

where the sum runs over the neighbours j of i, k is a, typical stiffness, i.e., 
the order of magnitude of the relevant Hessian components d^V/dridvj, and 
the ufs are the displacements with respect to an equilibrium configuration. 
Let us now introduce a continuous, coarse-grained displacement field u{r]t) 
and a typical interparticle distance oq, and substitute the former into Eq. 
in the overdamped limit 0, 

m du 
T dt 
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ka^V^u. 





FE 

MD 


4 / 


m 


(a) At = 2 



(b) At = 10 



(c) At = 1000 


Figure 3: Average displacement field induced by a shear transformation (at 
the centre of the cell), after a time lag At, for relatively low damping, C = 1 
(hence, r] = 0.726). The pink arrows represent the displacement vectors 
and the background colour indicates their norms. System size: (205cr^^)^, 
corresponding to 82^ finite elements. 

(Left) Finite Elements, het. iso.; {right) Molecular Dynamics. 
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FE 


MD 







(a) At = 2 



(c) At = 1000 


Figure 4: Average displacement field induced by a shear transformation, after 
a time lag At, for strong damping, C = 100 (hence, rj = 72.6). Refer to Fig. 
for the legend. 
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(a) C = 1 (»? = 0.726) 




(b) C = 10 ( 7 / = 7.26) 




(c) C = 100 ( 7 ? = 72.6) 

Figure 5: Average propagation radius as a function of time, for different 
damping magnitudes. 

{Red stars) MD data; {inverted cyan triangles) FE, het. iso.; {blue triangles) 
FE, het. aniso; {solid black line) FE, uniform system. 

{Left) log-log plot, {right) same data, in semi-logarithmic plot. 
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In this regime of negligible inertia, we thus obtain a diffusive equation for the 
particle displacements, consistently with the MD observations. 

5.2.2. Dissipative Particle Dynamics 

Very crudely, the DPD equations of motion (Eqs. [T]j^ are approximated 
by 


mil ~ ('^^(ujit) — Ui{t)) + k'^2{uj{t) — Ui{t)) 

U\i) {j\i) 

mu ~ + ka^V^u, (9) 

where ( = (oq). 

Equation]^ is a diffusion equation (on ii) only if the elastic force is neg¬ 
ligible, which will not be the case in practice. (More generally, Eq. can 
be solved with a space-time Eourier transform, or a joint Laplace-Eourier 
transform). 

It can also be seen in Eq. [^that, regardless of the value of C, the inertial 
term mil will always dominate at long enough wavelengths. In an unbounded 
system, this notably implies that the inertialess Brownian limit, which fea¬ 
tures an infinite transverse sound velocity, is singular. 


6. Effect of structural disorder in MD and in FE 

Let us now investigate the impact of elastic heterogeneity on the displace¬ 
ment field induced by an individual plastic event, he., the importance of 
fluctuations around the disorder-averaged response. 

The norm of the average displacement u (r; t) along a diagonal direction, 
at a long time lag At = 1000, is plotted in Fig. for C = 1 and C = 100, 
along with the associated standard deviation 5u, i.e., 


Su (r; t) = ^ (r; t) — u (r; t)]" 


where the brackets denote an average over the realisations of disorder. Inci¬ 
dentally, one may notice that, for ^ = 1 (Fig. [6a]), MD and FE do not coincide 
satisfactorily with respect to the average displacements, but this is mostly due 
to a loss of synchronization: the oscillations described in Section [5T] have not 
died out yet at this time lag and they are not exactly in phase in the dif¬ 
ferent systems. Had the true steady-state limit. At ^ oo, been reached (at 
the expense of much longer simulations), we would have expected much bet¬ 
ter agreement on the average displacements. This expectation is supported 
by the coincidence of the average displacements at At = 1000 under strong 


damping, for C. = 100 (see Fig. 6b), in which case dissipation is more efficient 


and the steady state is reached after fewer MD steps; indeed, in the linear 
regime probed here, the final state should be independent of the dynamics, 
hence of (. 

Regarding the fluctuations, the main result is that their order of mag¬ 
nitude is well reproduced by the FE simulations, both with isotropic blocks 
(het. iso., /ii = 112 ) and with anisotropic blocks (het. aniso.), although, quite 
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(a) C = 1. (b) C = 100. 

Figure 6: {Solid lines) mean value u and {dashed lines) standard deviation 
Su of the displacement norm along a diagonal axis Cdiag = ^ {^x + By), after 
a time lag At = 1000, as a function of the distance (in FE units). 

{Red) MD; {eyan) FE, het. iso.; {blue) FE, het. aniso. 


naturally, het. aniso. displays larger fluctuations than het. iso. Moreover, it 
is noteworthy that these corrections Su are roughly half as large as the mean 
reponse at a distance of, e.g., 50a aa- To avoid any misunderstanding on the 
possible nature of the fluctuations measured in MD, let us recall here that 
the centre of mass of the MD simulation cell is kept fixed, which prevents 
the variable global translations of the system that are sometimes observed 
otherwise (and which then dominate the fluctuations)!^ 

With regard to the spatial distribution of Su, colour maps of the relative 
fluctuations Su{r;t)/u{r;t) are presented in Fig. In regions with non- 
negligible displacements, i.e., u{r;t) ^ 10“^, the relative fluctuations are 
approximately homogeneous and tend to increase slightly with time. 

In conclusion to this section, taking into account the broad distribution of 
shear moduli in FE has enabled us to recover the fluctuations observed in MD. 
This further confirms the role of structural disorder on the redistribution of 
stress induced by a plastic event. In the last section, we go one step further by 
attempting to reproduce the individual, time-dependent response to a given 
plastic event in MD with the simple FE framework. 

7. Time-dependent response to a particular plastic event 

Even though the study of the propagation dynamics (Section and of 
disorder-induced fluctuations (Section]^ validates the FE method for (future) 
use in, e.g., mesoscopic rheological models, we would like to know whether 
the comparison can be pushed further. More precisely, can the FE routine 
describe the details of the elastic response in a particular conflguration? 


^When the centre of mass of the MD simulation cell is not kept fixed, the fluctuations 
5u measured in MD are significantly larger and their profile with respect to the distance r 
to the origin (dashed lines in Fig. is almost flat. 
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Figure 7: Colour map of the relative displacement norm fluctuations 
Su{r]t)/u{r]t) for ( = 1. The regions where u{r]t) < 10“^ are overlaid 
in light yellow. 
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To address this question, within the third type of FE mode, namely, het. 
aniso., the local shear moduli jii and /i 2 and the angle 6 of each FE macro¬ 
element [i.e., set of four adjacent elements) are directly extracted from the 
corresponding region in the MD system. Then, we compute the coarse-grained 
strain fielc|^ induced by shear transformations occurring at given position in 
the sample, an example of which is shown in Fig. 

Clearly, the MD response and its FE counterpart look alike and both ex¬ 
hibit the distinctive quadrupolar angular structure associated to the response 
in a uniform medium. However, are the disorder-induced fluctuations, i.e, 
the deviations from this average response, also similar in MD and FE? In an 
endeavour to answer this question, we have looked at the deviations in half a 
dozen particular conflgurations (not shown) and considered a couple of basic 
measures of similarity, but our results remain inconclusive in this respect: 
there is no quantiflable evidence that the disorder-induced fluctuations in a 
particular MD configuration are satisfactorily reproduced in FE. 


8. Conclusions 


In conclusion, we have extracted information about the local elastic con¬ 
stants of a binary Lennard-Jones mixture and the viscosity associated with a 


DPD damping scheme. Consistently with the findings of Mizuno et al. (2013), 


we have found that the local shear moduli are more broadly distributed (on 
a relative basis) than local bulk moduli. 

These elastic and viscous properties were used as input in a simple FE 
routine and an ideal shear transformation was artificially triggered in the (FE 
and MD) systems. 

We observed that the average time-dependent elastic response to this 
transformation in a disordered medium is similar to the propagation in a 
uniform medium and it is well reproduced in the FE simulations. However, 
fluctuations with respect to the average displacement held are considerable, 
with relative fluctuations of a few tens of percents. The approximate mag¬ 
nitude of these fluctuations is captured by FE simulations on heterogeneons, 
but locally isotropic systems. Refining the description by considering the 
elastic anisotropy on the mesoscale does not play a major role in this respect. 

It should however be stressed that, throughout our investigation, shear 
transformations were arbitrarily imposed, through an instantaneous displace¬ 
ment of particles (or FE nodes). However, in a bona fide simulation, the 
dynamics of shear transformations are determined by the system itself; two 
dynamical regimes can then be envisioned: 

(i) if inertia is negligible, the competitition between elasticity and viscosity 
sets the timescale of the rearrangement, r = rj/fi, 

(ii) if the rearrangement mostly consists in the damping of the inertial 
force (initially generated by elasticity), then the duration of a rearrangement 
is set by the inverse damping coefficient 


^In MD, local strains are computed after coarse-graining the displacement field on a 
grid similar to the FE one; note that the strain field is expected to be less sensitive to 
heterogeneities than the displacement field. 
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(b) At = 10 
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(c) At = 100 

Figure 8: Local strain field induced by a particular shear transformation at 
different lag times, for ^ = 1. 

{Left) FE, with an elastic configuration modelled on the MD system; {right) 
MD. 
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All in all, our method represents a powerful new framework for rheological 
models for amorphous solids, which improves on the traditional use of an 
analytical elastic propagator and the computation of the response by means 
of a Fast Fourier Transform, in that it accounts for structural disorder and 
inertial effects, whose impact has been underscored by jSalerno et al. (2012), it 
can be extended to arbitrary (in particular, confined) geometries, and it may 
include pre-existing local defects in the material, such as cracks. A further 
asset of this strategy is that, notwithstanding the enhanced capabilities of 
the algorithm, its complexity in terms of number of operations scales linearly 
with the number of blocks (or FE cells) for large systems, that is, with a 
scaling comparable to that of the Fast Fourier Transform routine. 
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Appendix A. Simplified Finite Element routine 

Bearing in mind our pursuit of minimalism, we choose a simple regu¬ 
lar square meshgrid, as sketched in Fig. [TJ If one assumes that the strain 
and stress fields are approximately uniform in each element, the following 
equations can be written between the (nodal) displacements {ux,Uy) and the 
(elemental) strains e, on the one hand, and the (nodal) forces (/|', /®') and 
the (elemental) stresses cr, on the other hand: 
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where the nodes of the element have been numbered from 0 to 3 connter- 

o o ( n'\ 

clockwise, starting from the bottom left corner, viz.^ and ni denotes 

the displacement along x at the (0) node, etc. Here, we have used condensed 
notations for the 2D strains and the stresses, viz., 


e = 


^XX ] 

i 

/ (7^^ 

I ^ XX 

^yy 

and a = 

^ty 

V^^xy ) 

f 

V \/2< 


and the matrix B is given by 
-1 


B = 1/2 


1 1 -1 
- 1-11 1 

-IA /2 -1/^ -I/V 2 -V^ 


Notice that our simplified FE method is close to a Finite Volume method, 
in practice. The prefactors have been introduced with foresight (see Sec¬ 


tion 3.2) and the “minus” sign preceding B in Eq. A.l should not come as a 
surprise if one recalls that is the force exerted by the element on node 




Contrary to traditional EE codes, the mesh will here remain static, he., 
not be distorted owing to the material deformation. 


Annexe A.l. Elastic force-displacement matrix 

The objective is now to rewrite Eq. |^in terms of nodal displacements and 
forces in order to arrive at Eq. 

To relate the nodal displacements and the nodal forces in each element, 
we make use of the constitutive equation of the material. 

To start with, the elastic contribution is governed by Hooke’s law, which 


reads, in condensed notations (Tsamados et ah, 2009), 


(T = C 


(A.2) 


where C is a 3 x 3 real matrix. Substituting from Eq. |A.1[ one obtains the 
local relation between the forces exerted on the nodes by the material element 
under consideration and the displacements at the nodes, viz., 
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(A.3) 


To proceed, the loeal elastic force-displacement matrices K = —B^CB 
are assembled into a global elastic force-displacement matrix K, viz., 
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where the bold superscripts refer to the global labels used in Fig. by oppo¬ 
sition with the elemental labels used in Eq. |A.3[ Here, /C is a sparse 2N x 2N 
matrix. 


Appendix A.2. Viscous force-velocity matrix 

The foregoing derivation relies on the linear relation connecting local 
strains and elastic stresses. Thus, it can straightforwardly be extended to 
the viscous stresses, insofar as they are linearly related with the local strain 
rates, viz., 

^diss _ ^diss . ^ 

Globally, the viscous force-velocity relation reads 


/ y*diss(N—1) \ 
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where the 2N x 2N matrix H has been assembled from elemental matrices 
of the form 


Appendix A.3. Inertial force-acceleration matrix 

Finally, we must express the inertial forces, that is to say, the matrix Ad 
in Eq. The convected part of the material derivative of the velocity, namely, 
V ■ (Vu) , which scales with for elements of unit size, is neglected. 

We compute the inertial forces directly at the nodes. In other words, each 
node is assigned a mass tuq = pVo, where Vq is the elemental volume (he., 
area). Accordingly, the lumped-mass matrix Ad connecting the accelerations 
at the nodes to the inertial forces at the nodes is a 2N x 2N matrix with mo 
on the diagonal, he., 

/ mo \ 

Ad = 

\ rno J 

Below, we detail the steps and approximations that bridge the gap between 
the Continuum Mechanics formulation of Eq. [^and the following EE problem, 
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elasticity 


viscosity 


inertial force 

where the Ux^’s and u^^’s are the displacements at the nodes i G {0,. 
of a regular mesh. 


,N-1} 
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Appendix A.4.- Discretisation of the dynamics 

A central difference scheme is nsed to discretise Eq. |A. 5 | in time, viz., 


5uiQ = 

2dt 

5uiQ = (Wi) + Suit ,) - 25u (4) ^ 


(A. 6 ) 


where tn-i, tn, and tn+i refer to consecutive time , separated by a fixed time 
step St. 

After insertion into Eq. A.5, provided that Su{tn-i) and Su{tn) are known, 
the displacements at the next time step 5u{tn+i) are straightforwardly ob¬ 
tained by inverting a matrix. The advantage of using a static meshgrid is 
that this matrix is then constant and, accordingly, can be inverted once and 
for all at the beginning of the simulation. 


Appendix A.5. Biperiodic boundary conditions 

We implement biperiodic boundary conditions by connecting the leftmost 
nodes of the system to the rightmost ones (see Eig. [^, and the top row to 
the bottom one. 


Appendix B. Relation between the intrinsic macroscopic viscosity 
and the microscopic damping coefficient 

In MD, the damping magnitude is set by the coefficient ( in the expression 
of the dissipative force (Eq. [^, whereas it is set by the viscosity 77 in EE. 
In order to match the damping in both simulations, we must connect the MD 
dissipative force to the viscous stress in EE, namely, = 2pe (see 

Eq.|5- 

To this end, we consider a pure shear situation, in which particles are 
strictly advected by the flow 

v(r) = e ■ r 

with e = (cy ® 63 , + 63 , ® Cy ). 


On the one hand, in MD, the microscopic dissipative stress on particle i 
(of volume Vq) is obtained with the help of the Irving-Kirkwood formula, viz., 

O-(n) = 

3 

= -C (Oi) ® Vij. 

j b 

Eocusing on the xy-component of the stress and setting Vi as the origin of the 
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frame, i.e., Ti = 0, for convenience, we get 
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Here, n is the average number density of the system and g(r) is the (alledgedly 
isotropic) pair correlation function. Equation B.l expresses the stress in a 
volume of space occupied by a particle; elsewhere the stress is zero. Therefore, 
the average stress in the material reads 


'xy 
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On the other hand, in FE, the shear stress simply obeys a^y = 2gexy. 
It immediately follows that 





g{r)w^ (r) r^dr. 


(B.2) 


If decreases fast (but smoothly) and the particles are hard and dense 


enough, so that g (r) exhibits a sharp peak at r = Uq, the viscosity in Eq. B.2 
can be further approximated as 
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where Zc is the coordination number, he., the number of first neighbours (at 
a distance r ~ uq). 

Equation |B.2| is valid for a one-component system, but the extension to 
binary mixtures, of components A and B, is straightforward; with transparent 
notations, the viscosity reads 
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In the considered Lennard-Jones system, this leads to 77 = 0.726 C- 


Appendix C. Determination of the local stiffness tensors 

With our condensed notations for the stress and strain tensors (Eq. [^, 
the macroscopic stiffness tensor of an isotropic material of bulk modulus K 
and shear modulus /r reads 

/ K + n K — n 0\ 

C = K-fi K + fi 0 

\ 0 0 2fi J 

In comparison, local stiffness tensors display rather unusual properties. 
To grasp the meaning of their (lack of) symmetries, some brief general con¬ 
siderations abont elasticity and deformation are in order. 

Suppose that a small macroscopic strain e is applied to a sample and focus 
on a mesoscopic region S. The local linear strain tensor e is defined as the 
symmetric tensor that best matches the displacements of the particles in S 
due to the applied strain. Only if the deformation is strictly affine over the 
whole sample do the local strain tensors equate to e. 

Because, for a given short-range interparticle potential, the local stress 
cr results from the local configuration of particles, it is reasonable (bnt not 
strictly necessary) to suppose the existence of a function / such that 

cr = f{e). 

Let us write the first-order Taylor expansion of /, provided that it exists. 


c^a /3 - 0-^0 = Ca/SysejS + O (||e||^) , (C.l) 

where a, /S E {x,y} and is the quenched stress in the original configura¬ 
tion. With condensed notations, Eq. C.l turns intc|^ 
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+ 0(||€f). (C.2) 


The affine strain-local stress approximation consists in replacing the compo¬ 
nents of e on the rhs of Eq. |C.2| with those of the affine strain e, in order to 


®As a minor technical detail, note that, because the tensorial multiplication 
involves a summation on both €xy and components of the second-rank tensor C 

may not exactly equate to their counterparts in the fourth-rank tensor for instance, 

Cxy,xy — ‘^^xyxy ’ 


27 










determine C more easily. For subregions of size larger than 5cr^^, |Mizuno| 


et al. (2013) showed that this approximation is quite reasonable, although it 


slightly underestimates the spatial fluctuations of the elastic constants. On 
the other hand, should the local stress on the Ihs be computed for a local 
deformation equal to e, i.e., should the system not be allowed to relax to the 
energy minimum after the application of the affine strain e, then we would 
obtain the so-called Born term C^, which largely overestimates the stiffness 
of the disordered material ( Mizuno et al.f 2013). 

For the time being, all components of the second-rank stiffness tensor C 
are independent. But, if the local stress derives from a (twice differentiable) 
local strain-energy density e, i.e., 

de 

^a(3 


de, 


then 






d^t 


dtajide^s 

It immediately follows that = C-jSa^'i this symmetry property is trans¬ 

ferred to the second-rank tensor C (thanks to the carefully chosen \/2 prefac¬ 
tors in Eq. C.2). Indeed, Tsamados et al.| ( |2009 ) observed numerically that, 
for coarse-graining regions larger than 5 Lennard-Jones particles in diameter, 
assuming a symmetric stiffness matrix C creates an error of less than 1% 
on the local stress evaluations. In the MD system under consideration, we 
quantify the asymmetry of the mesoscopic stiffness matrices, computed over 
regions of size a = Scr^^, with the following measure: 


AC 


V AC?, with AC = C - 
z—/ 2 

\ 

\ {xx,yy,xy} 


What should ||AC|| be compared with? At first sight, the answer would be 
||C||, but the latter is dominated by large symmetric terms involving the bulk 
modulus K ss 100. Thus, on second thoughts, it appears more informative 
to remove the terms involving K; ||AC|| should then be compared to, e.g., 
(Tr (C) — 2K) = 4 (/r), with (/i) = 18.8. From the histogram of ||AC|| values 
plotted in Fig. |C.9^ , it transpires that deviations from symmetry in C are not 
strictly negligible, but symmetry may nevertheless be a decent approximation. 


To further reduce the number of local parameters, the isotropic contrac¬ 
tion/dilation vector (^ 2 / 2 ^ 2 / 20 )^ is supposed to produce an isotropic com¬ 
pression and, thus, to be an eigenvector of C, ergo 


n — —C 

^xy^xx ^xy^yy 

CxX,XX etyy^yy 


The assumptions of tensorial symmetry and isotropic response to contraction 
come down to projecting C onto a matrix of the form 

/ a S /^ \ 

C'= I h a —j3 I with a, 5,/?, u G M, (C.3) 
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(a) Asymmetry. 


(b) Discrepancy with the matrix C' given in 
Eq.iu:! 


Figure C.9: Histograms of the approximation errors made when supposing 
that the local stiffness tensors C are (a) symmetric, (b) of the form given in 
Eq.[C3 


where a and P will be the averages of the pairs {Cxx,xxi Cyy^yy) and {Cxy,xxi —Cxy,yy), 

is 


lA'Cll = lie - C'l 


respectively. The approximation error, quantified by 
plotted in Fig. |C. 9 | 3 . As expected, the deviations are somewhat larger than 
were C only symmetrised, but they remain under control. 


For each matrix C', we compute the eigenvalues ci ^ C 2 ^ C 3 and define: 

- the small local shear modulus fii = Cil2, 

- the large local shear modulus /i2 = C2/2, 

- and the bulk modulus is iF = C 3 / 2 . 

The distributions of these local elastic constants are presented in Fig. [CTO 
and their mean values and standard deviations are summarised in Table [H It 
should be noted that the average eigenvalues of the projected tensor C' differ 
by 10% or less from the eigenvalues of the full local stiffness tensors C. 

The components of C' can then be rewritten as follows 


a = K + iJ,2 cos^ 26 + /ii sin^ 26 
5 = K — 112 cos^ 26 — Hi sin^ 26 
P = (1^2 - 

V = 2/^2 sin^ 26 + 2(11 cos^ 26 


where the angle 6 has been defined in Section 3.2 
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Ail /^2 K 


(a) 


(b) 


(c) K 


Figure C.IO: Histograms (number of counts) of the measured values of the 
local elastic constants /ri, [ 12 , and K in subregions of size ^uaa x in the 
MD system. 
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